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Formation of collapsing cores in subcritical magnetic clouds: 
three-dimensional MHD simulations with ambipolar diffusion 

Takahiro Kudoh^ and Shantanu Basu^ 
ABSTRACT 

We employ the three-dimensional magnetohydrodynamic simulation including 
ambipolar diffusion to study the gravitationally-driven fragmentation of subcrit- 
ical molecular clouds, in which the gravitational fragmentation is stabilized as 
long as magnetic fiux-freezing applies. The simulations show that the cores in an 
initially subcritical cloud generally develop gradually over an ambipolar diffusion 
time, which is about a few xlO^ years in a typical molecular cloud. On the other 
hand, the formation of collapsing cores in subcritical clouds is accelerated by su- 
personic nonlinear flows. Our parameter study demonstrates that core formation 
occurs faster as the strength of the initial flow speed in the cloud increases. We 
found that the core formation time is roughly proportional to the inverse of the 
square root of the enhanced density created by the supersonic nonlinear flows. 
The density dependence is similar to that derived in quasistatically contracting 
magnetically supported clouds, although the core formation conditions are cre- 
ated by the nonlinear flows in our simulations. We have also found that the 
accelerated formation time is not strongly dependent on the initial strength of 
the magnetic field if the cloud is highly subcritical. Our simulation shows that 
the core formation time in our model subcritical clouds is several xlO^ years, due 
to the presence of large-scale supersonic flows (~ 3 times sound speed). Once a 
collapsing core forms, the density, velocity, and magnetic field structure of the 
core does not strongly depend on the initial strength of the velocity fluctuation. 
The infall velocities of the cores are subsonic and the magnetic field lines show 
weak hourglass shapes. 
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Introduction 



Magnetic fields in molecular clouds play an important role in the early stage of star for- 
mation. Particularly, magnetic fields in molecular clouds can regulate the cloud collapse and 
fragmentation process. The important parameter for these processes is the mass-to-flux ratio 
M/$, where M is the mass and $ is the magnetic flux of the cloud. The mass-to-fiux ratio 
represents the relative s trength of gravity and the magnetic field. There exists a critical mass- 



to-flux ratio fM/^l^Ht (jMestel fc Spitzer Ill956l : IStrittmatter Ill966t iMouschovias fc Spitzer 
1976 : Tomisaka. et al. 19881 ). If M/$ > (M/$)r-vit. the cloud is supercritical and is prone to 
collapse. On the other hand, if M/ $ < (M/ $)crit5 a cloud is subcritical and cannot collapse as 
long as magnetic fiux-freezing applies. A similar condition M/<l> < (M/$)crit = l/(27rG'^/^) 
is required for stability against fragmentation of a n infinite uniform layer that is flattened 
along the direction of a background magnetic field (INakano fc Nakamura Ill978l ). where G is 
the gravitational constant. 



Mestel fc Spitzer I (Il956[ ) pointed out that even if clouds are magnetically supported, 
ambipolar diffusion will cause the support to be lost and collapse will begin. More generally, 
a subcritical cloud undergoes a gravitationally-driven fragmentation instabili ty that occurs 



on the ambipolar diffusion tirnescale rather than the dynamical timescale ( lL anger I Il978 



Zweibel I Il998l : ICiolek fc Basu 1 120061 ). The lengthscale of the instability is fundamentally 
the Jeans-scale in the limit of highly supercrit ical clouds, but can b e much larger when the 
mass-to- flux ratio is close to the critical value (jCiolek &: Basu II2006I ). The idea that the star 
formation is regul ated by anabipolar diffusion and the magnetic field h as been considered for 



many years (e.g. 



Shu et al. 1119871 . Il999l : IMouschovias fc Ciolek Ill999[ ) 



Magnetic field strength measurements through the Ze eman effect reve al that the mass- 
to-flux ratios of cloud cores are close to the critical value (jCrutcher II2004J ). These observa- 
tions are consistent with core formation driven by ambipolar diffusion in subcritical clouds. 
However, in order to assess whether the fragmentation is occurring in a subcritical or super- 
critical mol ecular cloud, rnagnet ic field measurements of the envelope or diffuse region are 
important. ICrutcher et al.l ( 120091 ) recently attempted to do this using the Zeeman effect, and 
argued that the mass-to-fiux ratio of the envelope (in four measured positions with beam 
diameters in the range of ~ 0.5 — 1 pc in the vicinity of four different cloud cores) is typically 
greater than that of the core. This cont r adict s the model of core formation in subcritical 
clouds. However, IMouschovias &: Tassis I (|2009[ ) contested their statistical analysis, and ar- 
gued that dropping the restrictive assumption that the magnetic field is constant across four 
measured envelope regions of differing morphology and density would result in the opposite 
conclusion, i.e., that the core mass-to-flux ratio is likely greater than that of the envelope. 
Future observational tests of this kind, but beneflting from increased amounts of source 
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i ntegration time and spatial resolution, can settle the current difficulties of interpretation 
( jCrutcher et al.ll2010[ ). The complementary method of measurements of polarized emission 
from dust grains, which reveal the magnetic field morphology in the cloud, generally show 
that the magnetic field in cloud cores is well ordered, and application of the Cha ndrasekhar- 



Fermi method yield 


s mass 


1998 




Girart et al. 


2006) 



on core scales (< Ipc) with thos e on large sca l es (> 200 pc) for several molecular clouds, and 



found a significant correlations. lAlves et al. I (120081 ) distinctly shows that the magnetic field 



is locally perpendicular to the large ffiamentary structure of the Pipe Nebula. These recent 
observations of polarized emission indicate that the magnetic field provides a dominant force, 
and that core formation may have been driven by ambipolar diffusion in subcritical clouds. 

Most nonlinear calculations of ambipolar-diffusion-driven evolution in subcritical clouds 
have focused on a single axisymmetric core, but some recent models focus on a fragmentation 
process that results in multiple cores. Nonlinear calculatio ns of ambipolar-diffusion - drive n 



fragmentation in subcritical clouds were ffist performed by llndebetouw fc Zweibel I (120001 ) 



( 


2004 


) and 


Basu et al. 


(2009a) 



ried out two-dimensional simulations of a magnetized sheet in the thin-disk approximation, 
which incorporates a finite disk half-thickness consistent with hydrostatic equilibrium and 
thereby includes the effect of magnetic pressure. They found that the fragment spa cing in the 



nonlinear phase agrees with the prediction of linear theory ( ICiolek fc Basu II2006I ). and that 
the the subcritical (or critical) model had subsonic infall, while the supercritical model had 
supersonic infall speed. The ffist three-dimensional simulation of the gravitational fragmen - 



tation with magnetic fields and ambipolar diffusion was performed by lKudoh et al. I (120071 ). 
which verified some of main results of the thin-disk models: the dichotomy between subsonic 
infall speeds in subcritical clouds and somewhat supersonic speeds in supercritical clouds, 
for example. 

Besides the magnetic field, the supersonic turbule nce in molecular clouds is also an im- 

port a nt component in the early stage of star formation ( MacLow &: Klessen 112004 : McKee &: Ostriker 
20071 ). The inclusion of superso nic turbulent initial con ditio ns to a fragmentation m odel in 
subcritical clouds was studied by lLi &: Nakamura I (120041 ) and lNakamura fc Li I (120051 ). adopt- 
ing the thin-disk approximation. They found that a mildly subcritical cloud can undergo 
locally rapid ambipolar diffusion and form multiple fragments because of the initial super- 
sonic motion in which the large scale wave mode dominates the power spectrum. The core 
formation occurs on the order of turbulei ice crossing time over the simulation box, which is 



comparable to the dynamical time scale. iBasu et al.l (l2009bl ) carried out a parameter study 
of the fragmentation regulated by gravity, magnetic fields, ambipolar diffusion, and nonlinear 
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turbulent flows. They confirmed the onset of runa way collapse in subcrit ical cloud is signif- 
icantly accelerated by the initial nonlinear flows. iKudoh &: Basu I (120081 ) also verified that 
the mode of turbulence-acceler ated magnetically - reKula ted star formation occurs in a fully 



three-dimensional simulation. iNakamura fc Li I (120081 ) also carried out three-dimensional 



MHD simulations of subcritical clouds including star star formation and bipolar outflows, 
and applied their model to the Taurus molecular cloud. Since three-dimensional simula- 
tions of subcritical clouds are resource-limited, large parameter studies have not yet been 
performed. 

In this paper, we study the fragmentation process in subcritical clouds, including am- 
bipolar diffusion, by fully th ree-d imensional simu l ations . Our study follows the previous 



ones of iKudoh et al. I (120071 ) and iKudoh &: Basu I (120081 ). We focus on the early stage of 
core formation and evolution, and do not follow evolution past the runaway collapse of a 
core. We carry out a parameter study by running a large number of models, and with higher 
spatial resolution than in our previous studies. We are especially interested in the effect of 
the large-amplitude nonlinear initial perturbations on the time evolution of the cloud frag- 
mentation, and discuss the mechanism of turbulence-accelerated star formation in subcritical 
clouds. 

Our paper is organized in the following manner. The numerical model is described in 
Section 2, results are given in Section 3, and a discussion of the results is given in Section 4. 
We summarize our results in Section 5. 



Numerical model 



The n umerical model used in t his paper is almost the same as that used by lKudoh et al 



( 120071 ) and iKudoh fc Basu I (120081 ). We solve the three-dimensional magnetohydrodynamic 
(MHD) equations including self-gravity and ambipolar diffusion. As an initial condition, we 
assume hydrostatic equilibrium of a self-gravitating cloud along the direction of an initially 
uniform magnetic field. We also assume an initially uniform structure perpendicular to the 
magnetic field lines. In this equilibrium sheet-like gas, we input a random velocity fluctuation 
at each grid point. 



2.1. Basic Equations 



We solve the three-dimensional MHD equations including self-gravity and ambipolar 
diffusion, assuming that neutrals are much more numerous than ions: 



— + i> ■ Vp = — pV ■ -u, 



(1) 



, 1„ 1 . „ „ , 

— + iv- W)v = — H J y< B - W, 

at p cp 



V X (i; X + V X 



cp 



-{j xB)xB 



j = — V X B, 

An 



V^ijj = AnGp, 



(2) 

(3) 
(4) 
(5) 



P 



(6) 



where p is the density of neutral gas, p is the pressure, v is the velocity, B is the magnetic 
field, j is the electric current density, ip is the self-gravitating potential, and is the sound 
speed. Here, we consider the molecular clouds whose number density is about 10^ — 10^ 
cm~^. In this case, the cooling time is generally much smaller than the dynamical time. 
Therefore, instead of solv ing a detailed energy equati on, we assume isothermality for each 
Lagrangian fluid particle (IKudoh fc Basu Il2003l 120061 ): 



dcs 
dt 



dcs 
dt 



+ V ■ Vc. = 0. 



(7) 



This means that each parcel of molecular cloud gas as well as surrounding warm gas (see 
Sec. 12. 2p retains its initial tempera ture. For the neutral-ion col lision time in equation ([3]) 
and associated quantities, we follow iBasu fc Mouschovias I (1l994l ). so that 



1.4 



m + rUn 

Pi{cTw)in 



where pi is the density of ions and (aw) in is the average collisional rate between ions of mass 
rrii and neutrals of mass m„. Here, we use typical values of HC0"'"-H2 collisions, for which 
(crw)in = 1.69 X 10~^ cm~^s~^ , m\lra^ = 14.4, and = 2 x 1.66 x 10~^^ g. We also assume 



that the ion density pi is determined by the approximate relation ( lElmegreen Ill979t iNakano 
1979h 
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where we assume k = 0.5 throughout this paper, and take i^' = 3 x 10 ^cm ^ as a typical 
value. By using equation IQ, equation ([8]) can be written as 

Tni = 7P"'/', (10) 

where 7 ~ 170.2 g^/^cm~^/^ s using the above typical values. Then, equation ([3]) becomes 



dB 

~dt 



V X (v X B) + V X 



1 



cp- 



3/2 



{3 xB)xB 



:iii 



2.2. Initial Conditions 



As an initial condition, we assume hydrostatic equilibrium of a self-gravitating one- 
dimensional cloud along the z-direction ( iKudoh &: Basu II2003I . l2006l ). The hydrostatic equi- 
librium is calculated from the equations 



dp 
az 



:i2) 



az 



(13) 



P 



(14) 



subject to the boundary conditions 

g,{z = 0) = 0, p{z = 0)=po, piz = 0) 



Po clo 



(15) 



where po and Cso are the initial density and sound speed at z = 0. If the initial sound speed 
(tempera ture) is uniform throughout the region, we have the following analytic solution ps 
found by ISpitzer I Jl942h : 

Psiz) = posech^iz/Ho), (16) 



where 



Ho 



CsO 



(17) 



is the scale height. However, an isothermal molecular cloud is usually surrounded by warm 
material, such as neutral hydrogen gas. Observations also show that the transition between 
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molecular gas and surrounding gas is quite sharp (IBlitz I Il99ll ). In order to simulate the 
situation, we assume the initial sound speed distribution to be 



where we take cL 



1 + tanh 



Zd 



(1^ 



2Hq, and = O.IHq throughout the paper. This equation 



gives a numerical model of the temperature transition of c^^/'^so with the transition 

length of Zd- By using this sound speed distribution, we can solve equations (fT2|) - (fT4l) 
numerically. The initial density distribution of the numerical solution shows that it is almost 
the same as Spitzer's solution for < z < Zc- 



We also assume that the initial magnetic field is uniform along the 2;-direction: 



Bz = Bq, = By = 0, 



(19) 



where Bn is constant. 



In this equilibrium sheet-like gas, we input a random velocity fluctuation at each grid 
point: 



= VaRm{x, y), Vy = VaRm{x, y) , 



0.0, 



(20) 



where Rm is a random number with spectrum vl oc fc" in Fourier space, where k = (kl+ky)^^'^, 
and the root mean square value is about 1. The -Rm's for each of Vx and Vy are independent 
realizations. In this paper, we apply the two distinct spectra of oc /c^ and oc The 
parameter Va shows the strength of the velocity fluctuation. Models with the same spectrum 
use the same pair of realizations of Rm for generating the initial perturbations. 



2.3. Numerical Parameters 

The constant G along with the two parameters Cgo and po that are associated with 
the initial state allow one to choose a set of three fundamental units for this problem. We 
choose these to be c^o, -f^o, and po for velocity, length, and density, respectively. These yield 
a time unit to = Ho/cso- The initial magnetic field strength introduces one dimensionless 
free parameter. 



Bi Bl 

the initial ratio of gas to magnetic pressure at z = 0. 



In the sheet-like equilibrium cloud with a vertical magnetic field, /3o is related to the 
mass-to-flux ratio for Spitzer's self-gravitating cloud. The mass-to-flux ratio normalized to 
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the critical value is 

/i5 = 2nG'/'^ (22) 

-Do 

where 

/oo 
psdz = 2poHo (23) 
-oo 

is the column density of Spitzer's self-gravitating cloud. Therefore, 

/3o = (24) 

Although the initial cloud we used is not exactly the same as the Spitzer cloud, (3q is a good 
indicator to whet her or not the magnetic fiel d can prevent gravitational instability in the 



flux- freezing case (INakano fc Nakamura Ill978l ). 



We also define a dimensionless parameter 

a = 7^2^, (25) 

which shows the strength of the ambipolar diffusion. We note that a can also be equated 
to Tni.o/to, where rni,o is the initial value of the neutral-ion collision time calculated from 
equation ( 1T0|) . By using the typical values in Section 2.1, we get a ^ 0.11, which is taken as 
a fiducial parameter in this paper. We also vary a as a free parameter in Section 3.3. 

Dimensional values of all quantities can be found through a choice of po and Cgo- For 
example, for Cso = 0.2 km s~^ and uq = po/rrin = 10^ cm~^, we get Hq ~ 0.05 pc, to — 
2.5 X 10^ yr, and Bq ~ 20 pG if f3o = 1. 



2.4. Numerical Technique 



I n order to solve t he equations numerically, we use the CIP method ( lYabe fc Aoki 

199ll : lYabe et al. Il200ll ) for equations ([T]), and ([7]), and the method of characteristics- 
constrained transport (MOCCT; Stone & Norman 1992) for equation ([3]), including an ex- 
plicit integration of the amb ipolar diffusion terrn . The combinatiori of the CIP and MOCCT 
methods is sum marized in iKudoh et al.l ( 1l999l ) and lOgata et al.l (120041 ). It includes the 
CCUP method (lYabe fc Wang I Il99ll ) for the calculation of gas pressure, in order to get 
more numerically s table results. The numerical code used for this paper is based on that of 
Ogata et all tooi ). 



In this paper, the ambipolar diffusion term is only included when the density is greater 



than a certain value, 



pc 



We let Per = 0.3po both for numerical convenience and due to 
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the physical idea that the low density region is affected by external ultraviolet radiation so 
that the ionization fraction becomes large, i.e., Tni becomes small (jCiolek &: Mouschovias 
19951 ). Under this assumption, the upper atmosphere of the sheet-like cloud is not affected 
by ambipolar diffusion. This simple assumption helps to avoid very small time-steps due to 
the low density region in order to maintain stability of the explicit numerical scheme. 

We use a mirror-symmetric boundary condition at ^ = and periodic boundaries in 
the X- and ^/-directions. At the upper boundary, aX z = z^^t = 4ifo; we also use a mirror- 
symmetric boundary except when we solve for the gravitational potential. This symmetric 
co ndition is just f or nu meri cal convenience. When the same boundary conditions were used 
by lKudoh et al. I ( 120071 ) and lKudoh &: Basu I (|2008[ ) for three-dimensio nal simulations, the re- 



sults were consiste nt with two-dimensional thin-disk simulations (e.g.. lBasu fc Ciolek 112004 



Basu et al.ll2009b( ). Hence, we believe that the boundary conditions do not affect the result 
significantly. The Poisson equation (|5]) is solved by the Green's function method to compute 
the gravitational kernels in the z-direct ion, along with a Fourier transform method in the x- 
and ?/-directions ( iMiyama et al.lll987bl ). This method of solving the Poisson equation allows 
us to find the gravitational potential of a vertically isolated cloud within \z\ < z^^^^. 

The computational region is \y\ < SttHq and < z < AHq. The number of grid 
points for each direction is {N^, Ny, N^, ) = (256,256,40). Sin ce the most unstable wavelength 
for no magnetic field is about AttHq ( IMiyama et al.lll987a( ). we have 64 grid points within 
this wavelength. We have also 10 grid points within the scale height of the initial cloud in 
the 2;-direction. The computational time of the fiducial model is about 70 hours of CPU 
time using four processors of the SX-9 in the National Astronomical Observatory of Japan. 
The maximum computational time, which occurs for the case of the highly subcritical cloud 
with small initial velocity fluctuations, is about 770 hours of CPU time. In the cases of the 
supercritical model, the computational times are about 30 minutes of CPU time. 



3. Results 

Table 1 summarizes the models and parameters for simulations presented in this paper. 
In the table, we have listed the values of the free parameters /3o, «, the form of the turbulent 
power spectrum, and Va- We also have listed the core formation time tcore, which is defined 
at the time when the density of the core reaches 100 po- Although the definition of tcore is 
determined from a practical constraint of the numerical simulation, in practice the center of 
the core is always supercritical at this time and the time evolution of the density shows the 
features of runaway collapse. In the models VO to V5, we change the amplitude of the initial 
velocity fluctuation Va, with the form of the turbulent power spectrum fixed at vl oc k~^. In 
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the models KO to K4, we change the amphtude of the initial velocity fluctuation, but with 
the form of the turbulent power spectrum fixed at oc . In the models AO to A4, we 
change the dimensionless ambipolar diffusion coefficient a for models with Va = 3.0cso or 
Va = O.lCso- In the models BO to B9, we change /So, the initial plasma beta at z = 0, which 
corresponds to the square of the initial mass-to-flux ratio, for models with Va = 3.0cso or 
Va = O.lCso- We are mostly interested in core formation in subcritical clouds, but also carry 
out some supercritical cases for comparison. 



3.1. General properties of a fiducial model 

We show the result of model V4 as a fiducial model, in which the initial turbulent velocity 
amplitude {va) is 3 times sound speed (cso), its spectrum is vf. oc k~^, the initial normalized 
mass-to-fiux ratio is about 0.5 (i.e., /3o = 0.25), and the dimensionless ambipolar diffusion 
coefficient has a typical value (a = 0.11). We suppose that these values are approximately 
typical in molecular clouds. 

Figure 1 shows the time snapshot of the logarithmic density colour map overlaid with 
the velocity vector field for model V4. The snapshot is at the end of our simulation, when 
the maximum density is greater than lOOpo- The upper panel shows the cross section at 
2; = 0, and the lower panel shows it at y = 20.7Hq. The maximum density is located at 
(x, y, z) = {—13.4:Hq, 20.7i/o) 0). A collapsing core is located in the vicinity of the maximum 
density. Figure 2 shows the time snapshot of the logarithmic plasma /3, the ratio of gas 
pressure to magnetic pressure, at the end of the simulation. It shows that f3 is greater than 
1 around the core. It means that once the collapsing core is formed, the mass-to-fiux ratio 
of the core is expected to be greater than 1 (/3 > 1), i.e., supercritical. The other collapsing 
cores are also formed in the region where plasma /3 is greater than 1 (the left-down region 
of the upper panel in Fig. 2). The left panel of Figure 3 shows the density, x- velocity, and 
plasma /3 along an a;-axis cut at y = 20.7ifo and z = 0. The right panel shows the density, 
2-velocity, and plasma /3 along a z-axis cut at x = —13AHq and y = 20.7ifo- In each panel, 
the velocity shows infall motion into the center of the core, though the x-velocity outside 
of the core (low density region) is affected by turbulent motions originating from the initial 
perturbation. The relative infall speed is subsonic and is about half of the sound speed at 
that time. The subsonic i nfall feature is consistent with previous stu dies of core formation 



in subcritical clouds (e.g.. ICiolek fc Basu II2OOOI : iBasu fc Ciolek Il2004j ). 



Figure 4 shows the time evolution of the maximum density at z = and /3 at the 
location of the maximum density. The density goes up to ~ 10 times greater than the initial 
density during the first compression of the cloud. Then, it rebounds and shows oscillations. 
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Eventually, the dense region goes into runaway collapse around t — 18to, which corresponds 
to about 5 X 10^ years. In the first compression, /3 also goes up to ~ 0.6. It gradually 
goes up and becomes greater than 1 when the runaway collapse occurs. As we discussed in 
the previous section, /3 roughly equals the square of the normalized mass-to-fiux ratio in the 
gravitational equilibrium state of gas with vertical magnetic field. The runaway collaping core 
is expected to be supercritical (/3 > 1). Figure 4, Figure 2, and Figure 3 demonstrate that 
^ is a good indicator of the mass-to- flux ratio, though the gas is not exactly in gravitational 
equilibrium at the last stage of our simulations. 

Figure 5 shows a close-up view of the core in the vicinity of (,x, y, z) = {~13AHo, 2{).lHo. 0). 
The iso-surface contour shows the logarithmic density and the lines represent the magnetic 
field. The core is located in a filamentary structure that was induced by the initial velocity 
fluctuation. The magnetic field hues show an hourglass shaped structure because of the infall 
motion into the center of the core. 



3.2. Effect of initial velocity fluctuations 

Figure 6 shows a time snapshot of the logarithmic density overlaid with velocity vectors 
for model VO. The model VO coresponds to the model for initial velocity amphtude {va) 
is 0.01 times sound speed (cgo), but other parameters are the same as the fiducial model 
V4. The snapshot is at the end of our simulation, when the maximum density is greater 
than lOOpo- When the initial perturbed velocity is much smaller than the sound speed, the 
core formation time {t 

core ^ 136^0) which corresponds to ~ 3 x 10^ years) is much longer 
than that (tcore ~ 18.7to) of model V4. The maximum density is located at {x,y,z) = 
{7.6Hq, ~11.9Hq,0). a collapsing core is formed in the vicinity of the maximum density. 
Figure 7 shows the time snapshot of the logarithmic plasma /3 at the end of the simulation 
for model VO. As we expected, the plasma /3 value is greater than 1 around the core ,which 
means the core is supercritical. In Figure 6 and Figure 7, core-like structures are gathered 
together on the side oi x > 0. This reflects the initial perturbation of A;"^ in which the larger 
scale contains greater energy in the velocity spectrum. The left panel of Figure 8 shows the 
density, x-velocity, and plasma /3 along an x-axis cut at y = 7.6Ho and z = 0. The right 
panel shows the density, 2;- velocity, and plasma /3 along a z-axis cut at x = —I.QHq and 
y = — ll.Oi^o- The velocity shows infall motion into the center of the core, though the center 
of the core shows a systematic speed in the positive x-direction of about 0.3 CgQ. The relative 
infall speed is also subsonic. Figure 9 shows the close up view of the core in the vicinity 
of {x,y,z) = (— 13.4ifo, 20.7i7o, 0). The iso-surface contour shows the logarithmic density 
and the lines represent the magnetic field. The core shows typical oblate-hke structure. The 
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magnetic field line shows a hourglass shaped structure and the deformation of the field line 
is similar to that of model V4 in Figure 5. 

Figure 10 shows the time snapshot of the logarithmic density overlaid with velocity 
vectors at the end of the simulation for model Kl. The model Kl coresponds to the model 
for the initial velocity spectrum is oc k^, but other parameters are the same as the fiducial 
model V4. The flat spectrum is possibly not consistent with observation, but it is used to 
compare with the result by the typical spectrum of oc When the initial perturbed 
velocity has a flat spectrum, the core formation time (tcore ~ 67to) is longer than that (tcore ~ 
18.7to) of model V4. The maximum density is located at y, z) = {—1.5Hq, —9.5Hq, 0) and 
a collapsing core is formed in the vicinity of the maximum density Figure 11 shows the time 
snapshot of the logarithmic plasma /3 at the end of the simulation for model Kl. We note 
that the /3 is greater than 1 around the core again. The left panel of Figure 12 shows the 
density, x-velocity, and plasma (3 along an x-axis cut at y = —Q.BHq and z — 0. The right 
panel shows the density, ^-velocity, and plasma (3 along a z-axis cut at x = —1.5Ho and 
y = —9.5Ho. The velocity shows infall motion into the center of the core, and the relative 
infall speed is also subsonic again. Figure 13 shows a close-up view of the core in the vicinity 
of {x,y,z) = {—1.5Hq,—9.5Ho,0). The magnetic fleld lines also show an hourglass shaped 
structure. Figure 3, Figure 8 and Figure 12 show the structure of collapsing cores are not 
so dependent on the initial velocity fluctuations, except that the outer region of the core in 
model V4 is strongly affected by the large-scale turbulent flow. 

The main difference of each model is the core formation time. Figure 14 shows the time 
evolution of the maximum density for models V4, VO and Kl. In the case of model V4, which 
has a A;~^ spectrum, the maximum density is strongly increased by the compression caused 
by the large scale supersonic ffow. In the case of model Kl, the compression occurs in the 
initial stage, but it is weaker than that of model V4 because the small scale perturbations 
contain more power than in model V4. In the case of model VO, with small- amplitude initial 
perturbations, the core formation time is about 7 times longer than that of model V4. 

Figure 15 shows the core formation time as a function of the strength of the initial 
velocity perturbation. As the perturbed velocity amplitude increases, the core formation 
time decreases for each type of initial spectrum. In the cases of the initial /c~^ spectrum 
(filled squares), the core formation times are systematically smaller than those of the initial 
A;° spectrum (open squares). In order to understand the physics of the core formation time, 
we plot the core formation time as a function of Ppeak in Figure 16, where Ppeak is defined 
as the value of the density peak during the first compression in the time evolution of the 
maximum density on z — 0. In the cases of the spectrum (filled squares), Ppeak is greater 
than in the cases of the spectrum (open squares), even when Va is relatively small. Figure 
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16 shows that the core formation time is shorter when the density peak Ppeak is greater, and 
indicates that it is nearly proportional to 1/ ^ ppeak- The density dependence is similar to 
that derived in quasistatically contracting magnetically s upported cores, assuming a fo rce 



balance of the magnetic force and gravity in the core (e.g., iMouschovias fc Ciolek Ill999l ). In 
our simulation, the density is increased by the compression caused by large-scale nonlinear 
flows, and a dense region undergoes a rebound and oscillations, which means the core is 
not in an exact force balance. Nevertheless, it is interesting that the core formation time is 
nearly proportional to 1/^Ppeafc- We discuss this subject further in Section 4. 



3.3. Dependence on ambipolar diffusion coefficient 

Figure 17 shows the time evolution of the maximum density for models V4, AO, Al and 
A2. The model AO, Al, and A2 correspond to the models for a = 0, a = 0.05, and a = 0.2, 
respectively, but other parameters of them are the same as those of V4. 

As a decreases, the core formation time increases. When a = (model AO), which 
means no ambipolar diffusion, a core was not formed during the calculation time, lasting 
until t ~ 268^0- Figure 18 shows the core formation time as a function of a. When the 
initial velocity perturbation is small, the core formation time shows nearly linear dependence 
on (open triangles). When the initial perturbation is nonlinear (filled triangles), the 
dependence on is a little steeper than a linear dependence. For model A2 {a = 0.2, 
Va = 3cso), the core formation occurs very rapidly, without showing an oscillation of the 
maximum density. When a is large enough, the supercritical region can appear quickly after 
the first compression, resulting in a significantly reduced core formation time. 

In figure 19, we show the time evolution of the total kinetic energies in the computational 
region for model V4, AO, Al, and A2. The oscillation of the energies comes from the strong 
compression and rebound of the clouds. The figure shows that the kinetic energies quickly 
decrease within a crossing time of the initial perturbed velocity across the computational 
region. The e-folding time of the energies depends slightly on a. As a increases, the e- 
folding time increases. In the case of a = (model AO), about 10 percent of the initial 
kinetic energy remains at the end of the calculation {t ~ 268to)- Basu & Dapp (2010) found 
that the significant portion (~ 1/2) of the initial kinetic energy remains in the flux- freezing 
case using a thin-disk approximation. Our result shows that the kinetic energy is more 
dissipative than in the thin-disk approximation. This is because of the different boundary 
condition in the vertical direction. We discuss this subject further in Section 4. 



- 14 - 



3.4. Dependence on initial mass-to-flux ratio 

Figure 20 shows the time evolution of the maximum densities for models V4, BO, Bl 
and B2, showing the difference between models with different /3o- The model BO, Bl, and 
B2 correspond to the models for /3 = 0.04, (3 = 0.09, and (3 = 0.36, respectively, but 
other parameters of them are the same as those of V4. As we have shown in Section 2, Pq 
approximately equals the square of the initial normalized mass-to-fiux ratio of the cloud. 
The figure shows that the core formation time is longer when /3q is smaller. The figure also 
shows that the initial density enhancement is smaller when /3o is smaller. Figure 21 shows 
the core formation time as a function of Pq. When is greater than 1, the initial cloud is 
supercritical. In the case of the supercritical clouds, the core formation occurs in dynamical 
time scale. When the initial velocity is large, a core forms during the first compression. The 
core formation time of the supercritical clouds does not strongly depend on the initial mass- 
to-fiux ratio (or Pq). When Pq is smaller than 1, the initial cloud is subcritical. The core 
formation times of a subcritical cloud are generally longer than those of a supercritical cloud. 
In the case of a highly subcritical cloud however, the core formation time does not strongly 
depend on the initial mass-to-fiux ratio. The transition occurs near the critical cloud, /3o ~ 1- 
This tendency is consiste nt with that of the g rowth time from the linear analysis using the 



thin-disk approximation (jCiolek fc Basu 1120061 ) and with that of the ambipolar diffusion tira e 



in quasistatically contracting magnetically supported cores (IMouschovias fc Ciolek I Il999l ). 
It is interesting that the same tendency is achieved even when the initial velocity fluctuation 
is nonlinear (filled circles). 



4. Discussion 



4.1. Timescale of core formation 



One of the problems for core formation in subcritical clouds is that the timescale of th e 



core formation is slower than t hat is expected from observations (e.g., Jiiina et al. I Il999l ) . 



In order to solve this problem, IZweibel I (120021 ) and iFatuzzo fc Adams I (120021 ) pointed out 
that the ambipolar diffusion rate can be enhanced i n a tu rbul ent medium with a fluct uating 
magnetic field. In a different way, iLi &: Nakamura I (120041 ) and lNakamura &: Li I (12005! ) found 
that the compression of the cloud by the l arge-scale turbul e nt flo w shortened the timescale of 
core formation even in subcritical clouds. iKudoh &: Basu I (120081 ) and this paper followed the 
latter idea using three-dimensional simulations without a thin-disk approximation. Here, we 
discuss how the compression shortens the timescale of core formation in subcritical clouds. 



Figure 16 indicates that the core formation time is nearly proportional to l/i/Ppeafc, 
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where Ppeak is the value of the density peak during the first compression in the time evo- 
lution of the maximum density at 2; = 0. The ambipolar diffusion time (Tad) is estimated 
from equation flTT]) as Tad ~ 47r(v27rG/«)p^/^L^/i?^, where L is the gradient length scale 
introduced by the turbulent compression. Because the compression by the nonlinear flow 
is nearly one- dimensional, the magnetic field scales roughly as 5 oc L^^ within the fiux- 
freezing approximation. The surface density E also scales as S oc L^^ . If the compression 
is rapid enough that ver tical hydrostatic e quilibrium cannot esta blished, then p oc S oc 
and Tad oc 

^5/2 

oc p ' (lElmegreen This means that diffusion 



can occur quickly if the turbulent compression creates small values of L or large values of 
p. This would lead to a rapidly rising value of /3 in Figure 4 at t/to ~ 1- If diffusion is so 
effective during the first turbulent compression that a dense region becomes supercritical, 
then it will evolve directly into collapse. Alternatively, the stored magnetic energy of the 
compressed (and still subcritical) region may lead to reexpansion of the dense region. If 
reexpansion of the initial compression does occur, then there is enough time for the vertical 
structure settle back to near-hydrostatic equilibrium. In this case, the density scales roughly 
lik e p (X T? (X L~^, and t ^ a oc L oc p~^/^. The ambipolar diffusion time is also estimated 
by iMouschovias k. Ciolek I (1l999l ) as Tad oc p^^/^ in quasistatically contracting magnetically 



supported cores, assuming a radial force balance between gravity and the magnetic Lorentz 
force. In our simulation, force balance is not exactly achieved, but the time average (of 
the oscillations) in the cores can be approximately in force balance. Since the first strong 
compression leads to rapid magnetic flux loss, a higher density region (than the initial value) 
is rapidly attained and then settles into an approximate force balance. Since the phase of 
continuing ambipolar diffusion in the oscillating high density region takes longer than the 
initial compression, the overall time scale of the core formation is approximately proportional 

— 1/2 

to Ppeak • Here, instead of the time-average density, Ppeak is used as a representitive density 
of the compressed gas, since it is complicated to determine the time-average density of the 
moving, oscillating and collapsing high density region. It is interesting that the relation 
attains even when the Ppeak is used. We expect that the average density of the compressed 
overdense region would nearly proportional to the peak density. 

The relation means that the density dependence of Tad oc p~^/^ is the same as that of 
the free fall time. Figure 16 shows that the actual time of the core formation is about 30 
times longer than the free fall time of gas with Ppeak when a = 0.11 and /3o = 0.5. The 
vale of the dimensionless coefficient (~ 30) is expected to be inversely proportinal to a 
(section 3.3), and not to be strongly dependent on the mass-to-flux ratio except when it is 
nearly critical (section 3.4). Even when the time-average density is used instead of the peak 
density, the value of the dimensionless coefficient would not be so different because of the 
weak dependence of p (oc p^^^^). 
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In some cases, collapse may occur during the first compression itself or soon thereafter, 

— 1/2 

and the approximate scaling of the core formation time with Pp^a^ will not hold. Whether 
or not this occurs depends not only on the strength of the turbulent compression, but also 
on the initial mass-to-flux ratio of the cloud (related to /9o) and the ambipolar diffusion 
coefficient (a). In model B2 (a subcritical model that is closest to critical, with normalized 
initial mass-to-flux ratio ~ 0.6), a reexpansion does occur but the collapse starts relatively 
quickly, in the third oscillation (Fig. 20). In model A2 {a = 0.2), which has the poorest 
neutral-ion coupling of all models, the collapse starts almost in the first compression (Fig. 
17). 

An overall conclusion is that the core formation occurs more rapidly than it would in the 
initial state due to the elevated value of p in the compressed but oscillating region. Since the 
core formation time is approximately proprtional to p~^^'^, 10 — 100 times desity enhancement 
is needed to get 3 — 10 times shorter core formation time th an that of the standard core 



formation model in subcritical clouds (e.g. iJijina et al. 1119991 ). The observed non-thermal 



velocity (3 — 10 times sound speed) is eligible to make such enhancemenr by the compression 
in the isothermal clouds, if its scale is larger than the Jeans length. 



4.2. Structure of cores 



Even though the core formation time is accelerated by the nonlinear flows, the density, 
velocity, and magnetic field structure of a core does not strongly depend on the initial 
strength of the velocity fluctuation (e.g.. Fig. 3, Fig. 8, and Fig. 12). The infall velocities 
of the cores are subsonic and the magnetic field lines show weak hourglass shapes. This 
result may be consistent with the result that the core formation time is proportional to 

—1/2 

Ppeak ■ Even when core formation is initiated by the initial turbulence, the core properties 
can be similar to the q u asista tically contracting magnetically supported cores discussed by 
Mouschovias fc Ciolek I (1999). The initial turbulence accelerates core formation, but it 
eventually dissipates in the dense r egion when the collapsing core is formed. Subsonic infall 
motions were found by lLee et al. I (120011 ) in an observational survey of starless cores. They 
found that the typical infall radii are 0.06 — 0.14 pc and that the infall speed lies in the range 
of 0.05 — 0.09 km s~^. These values are consistent with our results. The subsonic infall is 
the common featu re of core formation i n subcritical cloud s, which has been pointed out by 



ID ax isymmetric (jCiolek fc Basu Il2000l ) and 2D thin-disk (IBasu fc Ciolek II2004J : iBasu et al. 



2009a| ]bl) models. Because of the subcritical infall, the hourglass shap es are a little weake r 



than that of an initially supercritical cloud that shows trans-sonic infall (iKudoh et al. 1120071 ). 



More quantitative analysis of the field morphological difference and the infall speeds may 
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distinguish the mass-to-flux ratio in a molecular cloud. 



4.3. Energy dissipation 



In Figure 19, we showed the time evolution of the total turbulent kinetic energy in 
isothermal subcritical clouds. It shows that the kinetic energy dissi pates efficiently even i n 
the flux-freezing limit (a = 0). This is different from the result of iBasu &: Dapp I ( 120101 ). 
They found that the turbulent kinetic energy persists in an isothermal subcritical cloud 
in the flux- freezing limit. They showed that the energy dissipation occurs effectively only 
in the stage of the first nonlinear compression. After that, the turbulent energy persists, 
and oscillates about an average value. The average value is estimated to be about half of 
the initial turbulent kinetic energy for the same parameters as those of our fiducial model. 
On the other hand, our simulation shows that energy dissipation occurs effectively even 
in the flux- freezing limit; the energy dissipation occurs almost exponentially during several 
oscillations, although about 12 percent of the initial kinetic energy remains at the final stage. 
The dif ference is likely due t o the different vertical boundary conditions on the molecular 



cloud. iBasu &: Dapp I (120101 ) use a thin-disk approximation of the cloud, and adopt the 



current-free condition of the magnetic field outside of the disk, assuming the outside density 
is negligibly small. We study the problem by three-dimensional MHD and without a thin- 
disk approximation, but assume that the molecular cloud is surrounded by warm gas whose 
density is 10 times smaller than the molecular cloud. The computational region in our model 
is < 2; < 4:Ho, so does not cover a very large dynamic range of densities. Therefore, future 
three-dimensional simulations with a larger computational region along the 2;-direction and 
including lo w density gas outside o f the disk may make for more revealing comparisons with 
the result of lSasu fc Dapp I J201oh . 



4.4. Additional parameter surveys for numerical accuracy 

In addition to the parameter study in Table 1, we performed simulations of models 
with the same parameters as models V4 and VI, but different random realizations of the 
initial velocity fluctuations. The core formation time for different realizations varies by 
about 25% in these samples. The overall evolution and the result of the accelerated core 
formation by nonlinear flows were not altered by the different initial random realizations. We 
performed simulations with different spatial resolutions for the model V4. In the case of the 
highest resolution of {N^, Ny, N^) = (512,512,40), the core formation time is tcore = 14.2to- 
There is a slight tendency that the core formation time becomes a little shorter for the high 
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resolution cases (IKudoh &: Basu II2009I ). However, we should note that a realization of the 
random perturbations for the initial velocity fluctuation is also not the same for the different 
resolutions, because the random perturbations are input on all scales down to the grid scale. 
At the least, we can say that the result of the accelerated core formation by nonlinear flows is 
not altered significantly by spatial resolution. We also performed a simulation with a different 
boundary condition of magnetic fields for the model V4. Instead of the symmetrical boundaty 
at z = Zout, we used the vertical field condition for the magnetic field, i.e., = By = for 
z > ^^out- The overall evolution and the result of the accelerated core formation were not also 
altered by this. The core formation time for the different baounday varies by about 0.2% in 
this sample. 



5. Summary 

We have performed fully three-dimensional magnetohydrodynamic simulations of col- 
lapsing core formation in molecular clouds with subcritical mass-to-flux ratio, including 
ambipolar diffusion. Some of our major findings are as follows. 

• Core formation in subcritical clouds is generally slow. The core develops gradually 
over an ambipolar diffusion time. When the initial mass-to-flux ratio is 0.5 times a 
critical value, the formation time is about 3 x 10^ years for an initial midplane number 
density 10'* cm~^. 

• The core formation time is shortened by strong velocity fluctuations. When the average 
strength of the velocity fluctuation is 3 times the sound speed, the formation time is 
about 5 X 10^ years for the same cloud described above. 

• The core formation time scales as tcore oc 1/y/Ppeak, where Ppeak is the value of the 
density peak during the first compression in the time e volution of the maximum densi ty. 



The density dependence is similar to that derived by iMouschovias fc Ciolek I ( 1l999l ) . 



In the case of a highly subcritical cloud, the core formation time does not strongly 
depend on the initial mass-to-flux ratio even when there is strong velocity fluctua- 
tion. This tendency is c onsistent with the resu lt of the linear analysis of thin-disk 



approximation studied by ICiolek fc Basu I ( 120061 ) 



Once a core forms, the density, velocity, and magnetic field structure of the core do not 
strongly depend on the initial strength of the velocity fluctuation. The infall velocities 
are subsonic and the magnetic field lines show weak hourglass shapes. 
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Numerical simulations were done on the SX-9 at the Center for Computational Astro- 
physics in National Astronomical Observatory of Japan. S. B. was supported by a Discovery 
Grant from the Natural Sciences and Engineering Research Council of Canada. 
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Table 1. Model and Parameters 



A/TnHpl 


PO 




>i <3 r" 1" n m 




(•core / ''0 


1; (Ttn TYi n G 


vo 


0.25 


0.11 




0.01 


136 




VI 


0.25 


0.11 




0.1 


88.9 




V2 


0.25 


0.11 




1.0 


41.8 




V3 


0.25 


0.11 




2.0 


24.9 




V4 


0.25 


0.11 




3.0 


18.7 


fiducial model 


V5 


0.25 


0.11 




4.0 


16.1 




KO 


0.25 


0.11 


feO 


0.1 


141 




Kl 


0.25 


0.11 




3.0 


67.0 




K2 


0.25 


0.11 




5.0 


64.9 




K3 


0.25 


0.11 


k" 


8.0 


55.0 




K4 


0.25 


0.11 




10.0 


44.7 




AO 


0.25 


0.0 




3.0 


> 268 


no collapse 


Al 


0.25 


0.05 




3.0 


63.6 




A2 


0.25 


0.2 




3.0 


3.65 




A3 


0.25 


0.05 




0.1 


220 




A4 


0.25 


0.2 




0.1 


44.9 




BO 


0.04 


0.11 




3.0 


29.9 




Bl 


0.09 


0.11 




3.0 


24.8 




B2 


0.36 


0.11 




3.0 


6.94 




B3 


4.0 


0.11 




3.0 


1.01 


initially supercritical 


B4 


9.0 


0.11 




3.0 


0.97 


initially supercritical 


B5 


0.04 


0.11 




0.1 


114 




B6 


0.09 


0.11 




0.1 


108 




B7 


0.36 


0.11 




0.1 


68.6 




B8 


4.0 


0.11 




0.1 


6.93 


initially supercritical 


B9 


9.0 


0.11 




0.1 


5.99 


initially supercritical 


BIO 


oo 


0.11 




0.1 


5.38 


initially supercritical 



Note. — /3o is the initial plasma /3 at z = 0. a is a dimensionless ambipolar 
diffusion coefficient. Va is the amplitude of the initial velocity fluctuation, tcore is 
the time of collapsing core formation. 
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Fig. 1. — Logarithmic density contours at t = 18.7to for the model V4. The model V4 
is the fiducial model, in which the initial velocity amplitude Va = 3cso, its spectrum is 
v1 oc k~^, the initial normalized mass-to-fiux ratio is about 0.5 (i.e., /3o = 0.25), and the 
dimensionless ambipolar diffusion coefficient has a typical value {a = 0.11). The unit of time 
to is — 2.5 X 10^ yr, and Hq is ~ 0.05 pc. The top panel shows the x — y cross section at 
z = 0, and the bottom shows the x — z cross section at y = — 20.7iJo- Arrows show velocity 
vectors on each cross section. The unit of the velocity vector is three times the sound speed 



- 24 - 



Time = 18.7 
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X 

Fig. 2. — Logarithmic plasma (5 contours at t = 18.7to for the model V4. The top panel 
shows the x — y cross section at z = 0, and the bottom panel shows the x — z cross section 
at y = -20.7ifo- 
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Fig. 3. — Left: The density (solid line), velocity (dashed line) and plasma /3 (dotted line) 
along cut at y = 20.7 Hq and 2; = in the snapshot shown in Fig. 1 and Fig. 2. The 

x-positions are measured by offset from Xc = —IS.AHq, which is the maximum density point 
for the core. Right: The density, velocity and plasma (3 along a z-eods cut at x — —IS.iHo 
and y = 20.7ifo in the snapshot shown in Fig. 1 and Fig. 2. The hne styles are the same as 
those in the left panel. The unit of length Hq is ~ 0.05 pc. 
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Fig. 4. — Time evolution of the maximum density (solid line) at z = 0, and the evolution 
of plasma /5 at the location of maximum density (dashed line) for the model V4. The unit 
of time to is — 2.5 x 10^ yr. 




Fig. 5. — Isosurface of the logarithmic density, and the magnetic field lines near a core in 
the vicinity of ix,y,z) = {-13AHo,20.7Ho,0.0) in Fig. 1. 
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Fig. 6. — Logarithmic density contours at t = 136to for the model VO. The model VO 
coresponds to the model for initial velocity amplitude Va = 0.01 c^o, but other parameters 
are the same as the fiducial model V4. The top panel shows the x — y cross section at z = 0, 
and the bottom panel shows the x — z cross section at ?/ = — 11.9ifo- Arrows show velocity 
vectors on each cross section. The unit of the velocity vector is the sound speed c^o- 




Fig. 7. — Logarithmic plasma /3 contours at t = 136to for the model VO. 
shows the x — y cross section at ^ = 0, and the bottom panel shows the x — 
aXy = -11.9ifo- 
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Fig. 8. — Left: The density (solid line), a;- velocity (dashed line) and plasma (3 (dotted 
line) along cut at y = — 11.9i/o ^-nd 2; = in the snapshot shown in Fig. 6 and 

Fig. 7. The x-positions are measured by offset from Xc = 7.6Hq, which is the maximum 
density point for the core. Right: The density, 2;-velocity and plasma f3 along a 2;-axis cut 
at X — 7.6H0 and y — — 11.9i?o in the snapshot shown in Fig. 6 and Fig. 7. The line styles 
are the same as those in the left panel. The unit of length Hq is ~ 0.05 pc. 




Fig. 9. — Isosurface of the logarithmic density and magnetic field lines near the core in the 
vicinity of {x,y,z) = (7.6ifo, — 0.0) in Fig. 6. 
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Fig. 10. — Logarithmic density contours at t = 67.0to for the model Kl. The model Kl has 
an initial velocity spectrum vl oc k^, but other parameters are the same as in the fiducial 
model V4. The top panel shows the x — y cross section at 2 = 0, and the bottom panel 
shows the x — z cross section at y = — 9.5ifo- Arrows show velocity vectors on each cross 
section. The unit of the velocity vector is the sound speed Csq. 
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Fig. 11. — Logarithmic plasma /3 contours at t = 67.0to for the model Kl. The top panel 
shows the x — y cross section at z = 0, and the bottom panel shows the x — z cross section 
at y = — 9.5iJo- 
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Fig. 12. — Left: The density (solid line), a;- velocity (dashed line) and plasma /3 (dotted 
line) along cut at y = — 9.5i7o 2; = in the snapshot shown in Fig. 10 and 

Fig. 11. The x-positions are measured by offset from Xc = —l.bHo, which is the maximum 
density point for the core. Right: The density, 2;- velocity and plasma (3 along a 2;-axis cut at 
X — —I.5H0 and y — —9.5Ho in the snapshot shown in Fig. 10 and Fig. 11. The line styles 
are the same as those in the left panel. The unit of length Hq is ~ 0.05 pc. 




Fig. 13. — Isosurface of the logarithmic density and magnetic field lines near the core in the 
vicinity of {x,y,z) = {-1.5Ho, -9.5Ho,0.0) in Fig. 10. 
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Fig. 14. — Time evolution ol the maximum density a.t z = for the model V4 (solid line), 
VO (dashed hne), and Kl (dash-dotted line), showing the difference between models with 
different velocity perturbation. The model V4 is the fiducial model, in which Va — S Csq, its 
spectrum is oc the initial normalized mass-to-flux ratio is about 0.5 (i.e., — 0.25), 
and the dimensionless ambipolar diffusion coefficient has a typical value (a = 0.11). The 
model VO has Va = 0.01 Cgo, but other parameters are the same as those of model V4. The 
model Kl has a velocity spectrum vl oc k^, but other parameters are the same as those of 
model V4. The unit of time to is ~ 2.5 x 10^ yr. 
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Fig. 15. — Core formation time as a function of the initial velocity amplitude. The filled 
squares represent the results for the initial spectrum (model VO, VI, V2, V3, V4 and 
V5), and the open squares represent those for the initial spectrum (model KO, Kl, K2, 
K3 and K4). The unit of time is ~ 2.5 x 10^ yr. 
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Fig. 16. — Core formation time as a function of the density peak during the first compression 
in its time evolution. The filled squares represent the results for the initial k"^ spectrum 
(model VO, VI, V2, V3, V4 and V5), and the open squares represent those for the initial 
spectrum (model KO, Kl, K2, K3 and K4). The dashed line shows that the core formation 
time is nearly proportional to 1/ yJPpeak- The dotted line represents the free fall time of gas 
with density ppeak (= ^/ \/Gppeak) for comparison. The unit of time to is — 2.5 x 10^ yr. 
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Fig. 17. — Time evolution of the maximum density at 2; = for the model V4 (solid line), 
AO (dashed line), Al (dash-dotted line) and A2 (dotted hne), showing the difference between 
models with different a. The model V4 is the fiducial model, in which the dimensionless 
ambipolar diffusion coefficient has a typical value [a = 0.11). The model AO, Al, and 
A2 correspond to the models for a = 0, a = 0.05, and a = 0.2, respectively, but other 
parameters are the same as those of model V4. The unit of time to is ~ 2.5 x 10^ yr. 
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Fig. 18. — Core formation time as a function of a, the dimensionless ambipolar diffusion 
coefficient. Tfie filled triangles represent the results for an initial velocity fluctuation Va — 
3.0cs (model Al, A2 and V4), and the open triangles represent those for an initial velocity 
fluctuation Va — O.lCs (model A3, A4, and VI). The unit of time Iq is ~ 2.5 x 10^ yr. 
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Fig. 19. — Time evolution of the total kinetic energies for the models V4 (sohd hnc), AO 
(dashed hne), Al (dash-dotted line) and A2 (dotted hne), showing the difference between 
models with different a. The model V4 is the fiducial model, in which the dimensionless 
ambipolar diffusion coefficient has a typical value [a = 0.11). The models AO, Al, and A2 
correspond to a = 0, a = 0.05, and a = 0.2, respectively, but other parameters are the same 
as those of model V4. The unit of time to is ~ 2.5 x 10^ yr. 



-42 - 




Fig. 20. — Time evolution of the maximum density on 2; = for the models V4 (sohd line), 
BO (dashed hne), Bl (dash-dotted hne) and B2 (dotted hne), showing the difference between 
models with different /3o. The model V4 is the fiducial model, in which the initial normahzed 
mass-to-fiux ratio is about 0.5 (i.e., /3o = 0.25). The models BO, Bl, and B2 correspond to 
(3 = 0.04, 13 = 0.09, and f3 = 0.36, respectively, but other parameters are the same as those 
of model V4. The unit of time to is ~ 2.5 x 10^ yr. 
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Fig. 21. — Core formation time as a function of ^o, the initial plasma ^, aX z — 0. The filled 
circles represent the results for an initial velocity fluctuation Va — 3.0 (models BO, Bl, 
B2, B3, B4 and V4), and the open circles represent those for an initial velocity fluctuation 
Va = 0.1c, (models B5, B6, B7, B8, B9 and VI). The unit of time is ~ 2.5 x 10^ yr. 



